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Abstract 

A set of MapleV R.4/5 software routines for calculating the numerical evolution 
of dynamical systems and flexibly plotting the results is presented. The package 
consists of an initial condition generator (on which the user can impose quite general 
constraints), a numerical solving manager, plotting commands that allow the user to 
locate and focus in on regions of possible interest and, finally, a set of routines that 
calculate the fractal dimension of the boundaries of those regions. A special feature 
of the software routines presented here is an optional interface in C, permitting fast 
numerical integration using standard Runge-Kutta methods, or variations, for high 
precision numerical integration. 
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PROGRAM SUMMARY 

Title of the software package: Ndynamics. 
Catalogue number: (supplied by Elsevier) 

Software obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland (see applica- 
tion form in this issue) 

Licensing provisions: none 

Operating systems under which the program has been tested: Linux (RedHat 5.2 and Debian 2.0.34), 
Windows 95, Windows 98. 

Programming languages used: Maple V Release 4 and 5 and ANSI C. 
Memory required to execute with typical data: 32 Megabytes. 
No. of lines in distributed program, including On-Line Help, etc.: 1370. 
Keywords: Dynamical systems, fractal dimension, symbolic computing. 
Nature of mathematical problem 

Computation and plotting of numerical solutions of dynamical systems and the determination of the 
fractal dimension of the boundaries. 

Methods of solution 

The default method of integration is a 5th order Runge-Kutta scheme, but any method of integration 
present on the MAPLE system is available via an argument when calling the routine. A box counting 
method is used to calculate the fractal dimension of the boundaries. 

Restrictions concerning the complexity of the problem 

Besides the inherent restrictions of numerical integration methods, this first version of the package only 
deals with systems of first order differential equations. 

Typical running time 

This depends strongly on the dynamical system. With a Pentium II 450 PC with 128 Mb of RAM, 
the integration of one graph (among the thousands it is necessary to calculate to determine the fractal 
dimension) takes from a fraction of a second to several seconds. The time for plotting the graphs depends 
on the number of trajectories plotted. If there are a few thousand, this may take 20 to 30 seconds. 

Unusual features of the program 

This package provides user-friendly software tools for analyzing the character of a dynamical system, 
whether it displays chaotic behavior, etc. Options within the package allow the user to specify character- 
istics that separate the trajectories into families of curves. In conjunction with the facilities for altering 
the user's viewpoint, this provides a graphical interface for the speedy and easy identification of regions 
with interesting dynamics. An unusual characteristic of the package is its interface for performing the 
numerical integrations in C using a 5th order Runge-Kutta method. This potentially improves the speed 
of the numerical integration by some orders of magnitude and, in cases where it is necessary to calculate 
thousands of graphs in regions of difficult integration, this feature is very desirable. 
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LONG WRITE-UP 



1 Introduction 

A large number of problems in physics, chemistry, etc, can be represented in terms of 
a set of differential equations, which constitute a (very often nonlinear) dynamical sys- 
tem. Among these systems of interest, autonomous dynamical systems are particularly 
common. A ^-dimensional autonomous dynamical system can be represented by 

^F(X), (!) 

where X = (X 1 {t), X 2 (t), . . . , X d (t)) and F = (F X (X), F 2 (X), . . . , F d (X)). The X^t), 
with i = 1, . . .,d, represent field variables, usually related to quantities associated with 
the specific problem under consideration, t is a continuous parameter, which can often be 
considered as the time, while the functions Fj(X), % = 1, . . . , d, are general functions of 
their arguments. The <i-dimensional space constituted by the Xi(t) is called the system's 
phase space, with trajectories in this phase space representing the solutions of (0). It 
is possible to interpret the trajectories generated by all the initial conditions in phase 
space as analogous to the paths followed by the particles of a flowing fluid with velocity 
field X, where an overdot represents a derivative with respect to t. If F« = Fj(X,t), 
the system is said to be a non-autonomous dynamical system. However, introducing 
a new variable Xd+i = t, the (^-dimensional non-autonomous dynamical system can be 
transformed into a (d+l)-dimensional autonomous system. An autonomous system is said 
to be Hamiltonian if ([!]) represents the set of Hamilton's equations and the phase space 
has a simplectic structure. Otherwise we have a dissipative or non-Hamiltonian system. 
In this case, volumes in phase space are not conserved, or, equivalently, V-X = V-F ^ 0. 

Consider a general nonlinear autonomous dynamical system with dimension d > 3. 
Among these systems those exhibiting chaotic behavior have assumed particular interest 
and constitute a rule than an exception. Roughly speaking, chaos means extreme sen- 
sitivity to small changes in the initial conditions. Due to nonlinearity, fluctuations in 
the initial conditions of chaotic systems evolve such that they can completely alter the 
asymptotic outcome of the unperturbed trajectories in phase space. There are several 
ways of analyzing chaos in dynamical systems, such as the Liapunov exponents or the 
technique of Poincare surface of sections for Hamiltonian systems. (For good reviews 
and books on chaotic dynamical systems see, for instance |Lj, 00 We are particularly 
interested in certain features underlying chaotic dynamical systems, namely sets char- 
acterized by a non-integer dimension. Such sets are known as fractals. In the realm of 
dissipative systems, attractors with fractal properties have been called strange attractors, 
whereas in Hamiltonian systems there exist the strange repellers ||, identified as repellers 
with fractal properties. An important numerical task in investigating chaotic behavior in 
dynamical systems is therefore the determination of the fractal dimension 0. 

The presence of fractals is related to chaos, and so these structures are also associated 
with sensitivity to small changes in initial conditions. To show this relation, we follow 
closely the procedure of || . Let us consider a non-fractal basin boundary £ that separates 
two attractors A_ and A + , usually fixed points (this will be explained further in section 0). 
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An uncertainty in our initial condition X(to) = Xo that has strength e means that the 
initial condition lies somewhere in the region |X — Xo| < e. According to Figure 1, given 
an uncertainty e, the initial condition 1 always goes to the attractor A+. On the other 
hand, since the point 2 lies near the boundary E, due to the uncertainty e the actual 
trajectory may go to either attractor A + or attractor A_. The initial condition 1 is said 
to be e-certain, and the initial condition 2 is e-uncertain. Obviously, e-uncertain initial 
conditions are those which lie within a distance e of the basin boundary S. Let us call 
/(e) the probability of obtaining an e-uncertain initial condition, in the sense that it is 
the fraction of the d- dimensional volume of the phase space which lies within e of the 
boundary S. In the case of a simple non-fractal boundary, /(e) ~ e. This scaling law 
tells us that an improvement in the accuracy of the initial condition by, say, a factor 10 
(a reduction of e by 10), reduces /(e), and hence our probability of potential error, by a 
factor of 10. For the case of a fractal boundary the scaling law (see Appendix) 



is valid, where a is the uncertainty exponent. The box-counting dimension of the basin 
boundary, D , is given by @] 



with d denoting the dimension of the phase space. For a fractal boundary, D > d — 1, 
implying that a < 1, whereas for a non-fractal boundary, D = d — 1, and a — 1. Suppose 
that, for instance, a = 0.2 and /(e) ~ e a2 . Then, in order to reduce the probability of 
error /(e) by a factor of 10, a reduction of e by an order of 10 5 is needed. The improvement 
in prediction by increasing the accuracy of the initial conditions is seen to be less favorable 
the smaller a is. 

Related to the above discussion, this paper presents a set of software routines, im- 
plemented in MapleV R.4/5, together with an optional interface in C for fast numeri- 
cal integration, which, among other things, calculates the dimension of the fractal basin 
boundary. The main goal is to determine numerically the scaling law (0) and therefore the 
fractal dimension. Consider an initial condition X and perturb one of the coordinates, 
say Xk, in order to obtain two initial conditions X — e and X + e. By integrating nu- 
merically we can determine to which attractor or asymptotic state both of the perturbed 
initial conditions goes. We call the original initial condition uncertain if the perturbed 
initial conditions have distinct attractors. If a large number of initial conditions are cho- 
sen randomly in a given region of the phase space, it is possible to determine the fraction 
/(e) of these that are uncertain, where /(e) is the ratio between the number of uncertain 
conditions and the total number of initial conditions. This process can be repeated for 
several values of e. According to Bleher et. al. in j|, the quantities /(e) and /(e) are 
proportional, so that the parameter a can be determined by the scaling of /(e) with e. 

The package has the following features: 

• user-friendly software tools for fast numerical integration of both Hamiltonian and 
non-Hamiltonian dynamical systems, as well as a set of commands for the flexible 
presentation of the results; 

• an optional C interface for fast high-precision numerical integration; 




(2) 



D 



d — a 



(3) 
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• tools for determining sets of initial conditions that satisfy given constraints, ex- 
pressed, e.g., as g(X±, X 2 , . . . , Xd) > 0, where g(Xx, X 2 , . . . , Xd) is a given function 
of Xi, i = 1, . . . ,d; 

• the ability to perform combined symbolic and numerical studies by implementing 
these software tools in a symbolic computing environment. 

The paper is organized as follows: in section ^| we give a summary of the package's 
commands, followed by a detailed description^] of its most relevant commands, mainly 
Nsolve, for the generation of initial conditions and numerical integration of the system 
allowing 2D or 3D plots, View, for changing the ranges of all variables (including the time) 
and Fdimension, to determine the dimension of the fractal basin boundary; section ^] 
illustrates the application of various commands of the package to two examples with 
known results already presented in the literature; section £| analyzes the performance of 
the various integration schemes, and section [| contains some concluding remarks. 

2 The Ndynamics package 

Basically the Ndynamics package consists of a set of routines for calculating the numerical 
evolution of dynamical systems and plotting the results. The package generates initial 
conditions in a user-defined region, subject to constraints provided by the user. After 
calculating the evolution of each trajectory, the user may use its plotting commands to 
locate possible regions of interest — regions where the system presents chaotic behavior, 
and so on. Finally, a set of routines is available to calculate the fractal dimension of the 
boundaries using box-counting algorithms. 

Summary 

A brief review of the commands of the package is as follows:^ 

• Nsolve takes as arguments a dynamical system, the ranges of the initial conditions 
(from which it generates the initial conditions), the range for the time (independent 
variable), and the time step to be used when plotting the associated trajectories. 
The variables to be included in the final plot (which can be 2D or 3D) must also be 
defined; 

• View generates a plot from points already calculated, allows the user to modify the 
ranges for all the variables (including the time), and to swap between 2D and 3D 
plots; 

• Boxcount applies the box counting algorithm to determine the number of hypercubes 
of side e (see the introduction) covering the boundary; 

• Fdimension uses Boxcount for many values of e to calculate the fractal dimension 
of the boundary. 

2 Aside from this, the package itself contains on-line help in standard Maple format which can be 
viewed as the user's manual for all the routines. 

3 This subsection and the next one may contain some information already presented in the previous 
sections; this is necessary to produce a self-contained description of the package. 
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Description 



A complete description of the A T dynamics package's commands is found in the on-line 
help. Here we present the most relevant part of that description. 



2.1 Command name: Nsolve 



Feature: plot 2D/3D graphs calculated for a set of initial conditions, generated by the 
program, within ranges defined by the user. 

Calling sequence^: 

> Nsolve (sys, [dep_ranges, [indep_range ,plotstep] ] , Frame, 
optional_parameters) ; 



Parameters: 
sys 

dep_ranges 



indep_range 

plotstep 

Frame 



- a set of first order ordinary differential equations. 

- a list with the ranges for the dependent variables in the format: 
[x=a..b, y=c..d, ...], where x, y are the dependent variables and 
a..b, c.d, are the ranges in which the initial conditions will be 
generated. 

- the range for the independent variable: t=tl..t2, where t represents 
the independent variable. 

- a numerical value defining the interval (in the independent variable) 
for the plot. 

- a list with the variables to be plotted. 



Optional Parameters: 
initial 

random 

diagonal 

method=option 

method= [rk5C , number] 

Synopsis: 



- Forces the command to stop after generating the initial 
conditions. 

- Tells the program to generate the initial conditions 
randomly. 

- Makes the program look for the initial conditions 
along the main diagonal of the hypercuboid defined by 
the ranges of the dependent variables. 

- Indicates which numerical method of integration will 
be used, where option can be any one of the options 
available to dsolve in MAPLE (see the associated help). 

- Tells the program to generate the C code and manage the 
interface to the C routines, which use a 5th order Runge- 
Kutta method to integrate the equations more rapidly. The 
parameter number defines the integration step (constant). 



The Nsolve command is a part of the Ndynamics package, which, to brutally sum- 
marize things, is designed to calculate the fractal dimension of boundaries. The command 
that actually does this is Fdimension. In order to calculate this dimension, the process 



'In what follows, the input can be recognized by the Maple prompt >. 
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involves the demarcation of dynamically interesting regions of the system being studied. 
This demarcation is achieved by running Nsolve and looking for such regions by analyz- 
ing the resulting graphs. To facilitate this, Nsolve is designed to calculate and plot the 
trajectories associated with a given system of differential equations, subject to constraints 
and adjustments which can be tailored to the specific search at hand. In order to control 
and customize the output of Nsolve, we have introduced the global and environment 
variables described below. 

• Global Variables 

1. Coloring 

2. constraint 

3. initial_conditions 

4. number_ic 

5. initiaLtime 

The first item above controls in which color a trajectory will be displayed. Coloring 
is a boolean condition that is applied to the last calculated point of each trajectory 
so that different trajectories are displayed in one of two colors, corresponding to the 
cases where Coloring is true or false. Coloring defines the boundary, since it is 
defined to discriminate trajectories on either side of it. If Coloring is not assigned, 
then the command uses a default setting, which produces a plot with a certain color 
pattern. This is useful when the criteria to define the boundary have not yet been 
decided, but it is desirable, perhaps for the purposes of defining these criteria, to 
see the general flow of the trajectories. 

The second item, constraint, is a variable that defines the constraints imposed upon 
the initial conditions. All initial conditions generated will satisfy these constraints. 
This variable can be any valid boolean Maple condition. For instance, it may be two 
(or more) simultaneous conditions, such as 100 > x 2 — y 2 + z 4 > 20, where x, y, z are 
dependent variables of the differential system. This variable does not have a default 
value and can be left unassigned, in which case the command freely generates the 
initial conditions inside the hypervolume defined by the dependent variable ranges. 

Item number 3 is a facility that allows the user to supply the program with specific 
initial conditions. If this global variable is assigned, the program skips the part that 
generates the initial conditions and uses those supplied by the user. Its syntax is 
the same as that produced by the program with initial_conditions unassigned. 
We refer the reader to the on-line help that comes with the package for this syntax. 

In the case that initial_conditions is not assigned by the user, the global variable 
number_ic has to be defined, otherwise the program issues an error message. As 
the name suggests, it defines the number of initial conditions to be generated. 

Finally, the global variable initial_time allows the user to define the initial value 
for the independent variable to be anywhere in the interval stipulated for its rangef]. 

5 In the present version, when using the C interface, this facility is not available. 
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The integration of the trajectory is then performed in both directions, up to the 
limits of the time range. 

• Environment Variables 

To customize the interface between the program and the user, the following envi- 
ronment variables were introduced: 

1. _Env_print, 

2. _Env_store. 

_Env_print controls the level of printing during the execution of the command. 
When set to zero, no printing takes place; if _Env_print =1, (the default value) 
basic messages are printed; when set to 2, time-counting messages will be printed 
and, finally, for the value 3 all possible messages are printed. 

The variable _Env_store, controls in which file the results will be stored. The initial 
conditions, the range for the variables, the differential equation system, etc., are, 
by default, stored in a file called usedatal and the calculated points (trajectories) 
are stored (by default) in a file called usedata2. This arrangement is fine if only 
one process is running at a time. It will, however, result in file conflicts if the user 
wishes to perform simultaneous runs of the program from the same directory. The 
assignment of _Env_store by the user to the desired filenames avoids such problems. 

The arguments 

The first argument of Nsolve is a set containing the system of differential equations. In 
the present implementation, the program only deals with first order ordinary differential 
equations (ODEs). This is not a major drawback, since all systems of ODEs can be put 
into this form. The second argument is a list in which the first element is itself a list 
giving the ranges of the dependent variables. The second element is a list containing the 
range of the independent variable and the step (for that variable) for printing the graphs. 
The third argument gives the variables to be plotted (this can be two or three variables). 
Finally there are the optional arguments. 

The optional arguments 

The optional arguments can be given alone or in conjunction and in any order. By 
default, Nsolve takes care of everything for the user — it generates the initial conditions 
and fully integrates the system. If the optional argument initial is given then it stops 
after calculating (and saving in the appropriate file) the initial conditions. 

The optional arguments random and diagonal act similarly. There are thus three ways 
in which the command may generate the initial conditions in the hypercuboid defined by 
the dependent variable ranges. The default is to produce initial conditions homogeneously 
throughout the hypercuboid. However, if the argument random is used then the initial 
conditions will be randomly distributed inside hypercuboid. If, instead, diagonal is used, 
all conditions will be generated along the main diagonal f\ of the region. 

6 By "main diagonal" we mean the straight line inside the hypercuboid which joins the point defined 
by the minimum value of the ranges of the independent variables to the point defined by the maximum 
values of those ranges. 
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Many different algorithms can be used to numerically integrate the system. Maple 
provides a great variety of these and all of them can be used with our command, via 
the option method. For example, if method= classical [rk2], the second order Runge-Kutta 
method will be used []. 

Apart from the numerical methods available in the MAPLE system, our program has 
a further possibility: a C-based numerical integrator. The only additional requirement for 
this is that the machine where the program is being run has a C compiler. In this case, 
the integrator is a fixed-step fifth order Runge-Kutta scheme based on the procedure RKQC 
in ||. The idea behind this option is to gain speed over Maple's interpreted procedures, 
thereby combining the flexibility of the symbolic environment with the efficiency of the 
compiled C routine (see section |4] for comparisons). It is our intention to introduce more 
compilable algorithms in later releases. 

2.2 Command name: View 

Feature: This commands allows easy and efficient viewing of different parts of the tra- 
jectories calculated by Nsolve, allowing the user to zoom in and out, and to change the 
variables being viewed (if necessary the user can even swap between 2D and 3D plots). 

Calling sequence: 

> View(list_range) ; 

Parameters: 

list_range - a list with the new ranges of those (two or three) variables 

the user wants to view in detail. The format is: [x=a..b, y=c..d, ...], 
where x, y are (dependent or independent) variables and a..b, c.d, 
their ranges. 

Synopsis: 

The importance of this command lies in the search for the relevant regions for studying 
the system's behavior. The user generally starts by running Nsolve in a broad region of 
the variable space and then studying more closely the regions, for example, where chaos 
is suspected. The graphs that View produces (2D or 3D) do not need to display the same 
variables that Nsolve uses in its argument Frame. 

2.3 Command name: Boxcount 

Feature: this command analyzes the effect of a small perturbation on the initial conditions 
on the fate of a bundle of trajectories. 

Calling sequence: 

> Boxcount (epsilon, final_time, optional_parameters) ; 

7 Details on the possible numerical solving methods available (and their syntax) are available via the 
MAPLE help page for the command dsolve [numeric] . 
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Parameters: 
epsilon 
final time 



- defines the size of the perturbation around each initial condition. 

- indicates the final value of the independent variable. 



Optional Parameters: 
method=option 

method= [rk5C , number] 



- Indicates which numerical method of integration will 
be used, where option can be any one of the options 
available to dsolve in MAPLE (see the associated help). 

- Tells the program to generate the C code and manage the 
interface to the C routines, which use a 5th order Runge- 
Kutta method to integrate the equations more rapidly. The 
parameter number defines the integration step (constant). 

Synopsis: 

Suppose we have generated a set of initial conditions (see |2TT 
each of these initial conditions within a neighborhood of radius = e. 
close to the boundary (that separates the different end-point behaviors) the two perturbed 
trajectories will have different fates. Boxcount then counts the number of initial conditions 
for which the pairs of perturbed trajectories evolve to different attractors. 



i. Boxcount perturbs 
For points sufficiently 



2.4 Command name: Fdimension 



epsilon_range 
TF 

Nepsilon 



savef ile=f ilename 



Feature: Fundamentally applies Boxcount a user-defined number of times, with a different 
value of e each time. 

Calling sequence: 

> Fdimension (epsilon_range , TF, Nepsilon, optional_parameters) ; 

Parameters: 

defines the range of variation of e. 
indicates the final value of the independent variable, 
indicates how many values of e will be considered in 
epsilon_range. 
Optional Parameters: 

- Indicates in which file the results will be stored. The 
default value of filename is XY-f ile. Since 
Fdimension makes many calls to Boxcount, only 
the pairs of data e, and the logarithm of the 
ratio of the number of initial conditions on the boundary 
to the total number of initial conditions (see appendix |A]) 
are saved to the file. 

- Indicates which numerical method of integration will 
be used, where option can be any one of the options 
available to dsolve in MAPLE (see the associated help). 

- Tells the program to generate the C code and manage the 
interface to the C routines, which use a 5th order Runge- 
Kutta method to integrate the equations more rapidly. The 
parameter number defines the integration step (constant). 



method=option 



method= [rk5C , number] 
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Synopsis: 

Fdimension essentially implements the box-counting process in our program. It uses 
Boxcount to evaluate the number of points on the boundary for many values of e thus 
calculating the fractal dimension of the boundary (using the method described in ap- 
pendix [A]) . 

3 Examples 

To illustrate the routines presented in the previous section we consider the well-known 
Lorenz system |J: 

x = a(y — x), 

y = —y — xz + Rx, (4) 
z = xy — bz, 

where a, R and b are constant parameters of the problem. R plays a crucial role as far as 
the dynamics of phase space is concerned. A simple analysis shows the presence of three 
critical points (the points satisfying the conditions x = y = z = 0): the origin Po(0, 0, 0) 
and, provided R > 1, two further points are symmetrical with respect to the z-axis, 
P± = (±^Jb(R — 1), ±^Jb(R — 1), R — 1). The stability properties of these critical points 
depend only on R, so that we set o = 10 and 6 = 8/3. We summarize the behavior of 
the solutions of the Lorenz system for distinct values of R > 1 that will be important in 
the applications of the routines and commands we have described. For 1 < R < 24.74 the 
origin P is unstable and P± are attracting stationary solutions, therefore representing 
two possible asymptotic configurations (except for the set of trajectories of measure zero 
that stay in the neighborhood of Po). However, around the critical value, R c ~ 13.926, 
Po develops into a homoclinic point such that, above this value, the basins of attraction 
around P_ and P + are no longer distinct. This means that trajectories can cross backwards 
and forwards between the two before settling down. Another transition occurs when 
R w 24.74, for which P_ and P + become unstable. Lorenz considered R = 28 and 
obtained a remarkable behavior displaying for the first time a strange attractor. As we 
shall show, the dynamics are highly erratic: a given trajectory can spiral around one of 
the critical points, P + or P_, for some arbitrary period then jump to the neighborhood 
of the other critical point, spiral around that for a while and then jump back to the first 
one, and so on. 

In Figure 1 we show the form of the trajectories for the Lorenz system, with coefficients 
a = 10, b = 8/3 and R = 20.06. We are therefore inside the parameter range for erratic 
behavior discussed above. We see that, given enough time (in our case >50 units of 
time), the trajectories settle down around one of the two possible asymptotic points, P±. 
Figure 1 was obtained with 10 initial conditions, (xq, yo, zo), that were randomly chosen in 
the volume of phase space defined by 0.9 < xq < 1.1, 0.9 < yo < 1.1 and 21.9 < zo < 22.1. 
This is our first illustration of the application of Nsolve. Figure 1 was obtained using 
the set of commands below. First, we set the accuracy, read the program and define the 
system: 
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> Digits := 16; 



> with(Ndynamics) ; 

[Boxcount, Fdimension, Nsolve, View] 

> sigma:=10.0: b:=8.0/3: R:=20.06: 

> Lorenz : =-[dif f (x(t) , t)=sigma* (y (t)-x(t) ) , 

diff (y(t) ,t)=-x(t)*z(t)+ R*x(t)-y(t) , 
diff (z(t) ,t)=x(t)*y(t)-b*z(t)}; 

Note that the integers in (Q) have been swapped for floating-point numbers. This is 
important in the case of the parameter b to avoid truncated integer division occurring in 
the C subroutines generated. 

We then define the input: 

> x.range := 0.9..1.1: y.range := 0.9..1.1: z.range := 21. 9.. 22.1: 

> t.range := 0. .50: 

> plot_spacing := 0.03: 

> Frame := [x(t) ,y(t) ,z(t)] : 

> number_ic := 10: 

> time.step := 0.001: 

Now we can call Nsolve: 

> graph := Nsolve (Lorenz, [[x=x_range,y=y_range,z=z_range] , 

[t=t_range ,plot_spacing] ] , Frame, method=[rk5C,time_step] ) : 

The use of Nsolve above illustrates the application of the global variable Coloring. 
Here we have based the choice of color by specifying the approximate location (x coor- 
dinate) of one of the asymptotic points P±. Trajectories ending around this point are 
colored green, while trajectories ending around the other critical point are colored black. 
Using Display (graph) gives the result shown in Figure 1. 

In this paper we take as a measure of the degree of chaos of a dynamical system, such 
as the Lorenz system above, the fractal dimension associated with the possible different 
exit modes under small changes of initial conditions, as discussed in section |l]. A note 
here is in order about our criterion for choosing the intervals between the values of e (the 
variations for the initial conditions as explained in section |2.4|) . Given ranges for the 
initial conditions which define a hypervolume V- 1C in d dimensions, and distributing N- 1C 
initial conditions throughout this volume, the largest e cannot exceed the average distance 
I between the initial conditions, since otherwise the "perturbation" e for different initial 
conditions will overlap leading to loss of predictability of the fractal dimension. We take 
for / the minimum distance between two initial conditions for a uniform distribution, 
which is easily deduced to be 
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In all the numerical tests performed the largest value for e is of the order of I. The 
lowest value for e used is determined by the statistical error associated with the number 
of initial conditions N falling inside the fractal boundary. We then have that A 1 / 2 will 
give the statistical uncertainty associated with the probability / = N/N[ C of a given initial 
condition falling within the boundary J3J. The numerical tests performed have values of 
e and N ic such that A>400, thus assuring an estimated error (VN /N) of at most 5% in 
the determination of /. 

The fractal dimension can be readily obtained by using the command Fdimension. 
We begin by generating the initial conditions as explained in section 

> number.ic := 10000: 

> Nsolve(Lorenz, [[x=x_range,y=y_range,z=z_range] , 

[t=t_range,plotting_spacing] ] , Frame , initial) : 



Note that 10000 initial conditions are used. By omitting the optional parameter random 



(see section [2.1|) , these number _ic initial conditions are uniformly distributed throughout 
the interval volume V- lc = (0.2) 3 defined by the dependent variable ranges. The use of 
random initial conditions leads to no appreciable change in the results, which always 
remain within the statistical error defined above. Values of e ranging from 10 -5 to 10~ 4 
were used in Fdimension when computing the fractal dimension. An easy way to check 
the number of initial conditions, N, which are generated inside the basin boundary is by 
using Boxcount. For example, having generated 10000 initial conditions, we can use: 

> Boxcount (0.00001, 60, method= [rk5C,0 . 01] ) ; 

From the, 10000, points (that were testable), 2202, were close to the boundary. 

[2202,10000] 

tells us that, for the N ic = 10000 initial conditions and for e = 10 -5 , N = 2202 initial 
conditions are around the boundary, with a resulting statistical error of around 2%. If we 
decrease e, say to 10 -10 , we get instead N = 551, and the statistical error is of the order of 
N~ l l 2 th 4%, close to the 5% limit of acceptability which we have imposed. Accordingly, 
for smaller values of e we must increase the number of initial conditions in order to assure 
that we have at least the minimum number of initial conditions on the boundary necesary 
for statistical acceptability. With this done, the fractal dimension is then obtained by 
issuing the command: 

> Fdimension(0. 00001. .0.0001, 60, 10, method=[rk5C,0.01] ) ; 

where we are evolving the trajectories through 60 units of time and we are taking 10 
equally spaced points in the range ln(10 -5 ) to ln(10 -4 ). The results obtained are shown 
in Figure 2 in terms of the logarithm of the number of uncertain initial conditions, (those 
that change their fate, or "color" on the graph), against ln(e). The slope of the best-fit 
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straight line through the points gives the number a (see Appendix). The fractal dimension 
is then defined from (|3]) as Dq = d — a. We obtain for the fractal dimension the result 
Dq ~ 2.87 ±0.01, in complete agreement |7j with the known result for the Lorenz system 
for the choice of parameters o = 10, b = 8/3 and R = 20.06. 

In Figure 3 we once again compute the fractal dimension with Fdimension, but with 
10" 8 < e < 10" 7 . The result is the same up to statistical error. 

As a further check of our set of routines, we evaluate the fractal dimension for the set 
of parameters o = 10, b = 8/3 and R = 10, in the parameter region prior to the onset 
of transient chaos ||. Using Fdimension, with 5 x 10~ 3 < e < 9 x 10~ 3 (an interval 
that ensures the minimum statistical error), we obtain the result shown in Figure 4, from 
which we conclude that Do ~ 1.98 ± 0.05. Further refinements of this result for a larger 
number of initial conditions leads us to conclude that D = 2 and a = 1 and that, as 
expected, there is no chaos for this choice of parameters. 

As a demonstration of the use of View, we illustrate the basin-boundary in the two 
cases analyzed above. Figure 5 shows the case R = 10.00 while Figure 6 illustrates the 
results for R = 20.06. In both cases number _ic is taken to be 100. 

The arguments to View are [x = — 0.1. .0.1, y = —0.1. .0.1, z = 0.9. .1.1] (for Figure 5) 
and [x = 0.9. .1.1, y = 0.9. .1.1, z = 21.9. .22.1] (for Figure 6). From Figure 5 we easily 
identify a well-defined (non-fractal) boundary separating the sets of initial conditions 
whose associated trajectories approach P + and P_. In contrast we see in Figure 6 that 
neighboring initial conditions may lead to very different end-states. This is the fractal 
basin boundary discussed in the introduction. 

4 Performance 

To point out the advantages of our hybrid symbolic/numeric approach, which uses Maple 
to manage the source code generation and compilation for the number crunching, while 
maintaining Maple's flexibility, table [I] presents a comparison of the elapsed time taken to 
perform the same calculation using the C interface and performing the entire calculation 
in Maple, using some of its numerical integration routines. To obtain a fair comparison 
particular care was taken to use methods of the same order and adjusting parameters to 
produce solutions with the same precision. The following procedure was used: for a given 
set of initial conditions, namely x = 0.679149319354506, y = -0.5692394267519960 
, zq = 22.00017552807044), Maple's inbuilt high-precision Taylor Series intergrator was 
used to integrate system (|) from £ = 0to£ = ll until converegence was obtained to 
an accuracy of 13 decimal places; subsequently Maple's global variable Digits was set 
equal to 13 (higher values of Digits seemed to cause problems with Maple's integrator), 
and Maple's inbuilt 5th order Runge-Kutta-Fehlberg integrator, rkf 45, was used up to 
t = 11. This integration was found to be accurate to 4 decimal places. The step size in 
the C routine was then adjusted to give the same| precision, resulting in a step size of 
0.002. An average was then taken for 200 integrations. 

The results below were obtained with Maple V.4 running in Windows 98 on an AMD- 
K6 266 with 64Mb of SDRAM. The C compiler used was Delorie's implementation of 
GNU's gcc, available form www.delorie.com, with level 3 optimization. 

8 in fact, the C routine was slightly more precise 
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Integration 
Method 


C-interface 
rk5 


Maple 
rkf45 


Seconds per Trajectory 


0.1 


7.0 


Ratio to the Fastest Case 


1 


70 



Table 1: Comparative performances of Maple's inbuilt numerical integrators and our C 
interface for the Lorenz system. 

5 Conclusions 

We have presented a set of software tools which allow great flexibility for the analysis of 
dynamical systems composed of first-order ordinary differential equations, Even though 
the construction of these programs was mainly motivated by the computation of fractal 
dimensions of basin boundaries of chaotic systems, this is not their only use, as was 
highlighted in the previous section. Large regions of phase space can be readily studied 
with a fairly large number of initial conditions, allowing the user to quickly find regions 
of interesting behavior. 

The optional numerical interface allows for faster numerical integration than with 
Maple's internal routines. The ease of use of the Maple software is thus combined with 
a faster numerical interface (here implemented in C, though this can be easily modified 
to use routines written in other languages such as Fortran), thereby allowing intensive 
numerical study of dynamical systems. 

A Fractals 

A.l Fractal Dimension 

The idea of a fractal is hard to define precisely. But, when we refer to a set F as a fractal, 
we typically have the following in mind: 

1. Fhas a "fine structure", i.e., detail on arbitrarily small scales; 

2. F is too irregular to be described in traditional geometrical language, both locally 
and globally; 

3. often F has some form of self-similarity, perhaps approximate or statistical; 

4. usually the "fractal dimension" of F (defined in some way) is greater than its "topo- 
logical dimension"; 

5. in most cases of interest, F is defined in a very simple way, perhaps recursively. 

There are many ways to define the dimension of a fractal. An important one is 
the Hausdorff dimension [|J - one possible generalization of the "primitive" notion of 
dimension. 

Suppose that we have a hypercube of side a. Its hypervolume is 
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H v = a d , (6) 

where d is the dimension of the hyperspace. 

If we divide the hypervolume into N hypercubic cells of side e, we have that 

H v = Ne d . (7) 

Dividing (|7]) by @ we obtain 

l = N ( e - 

\a 

Defining 5 = e/a and noticing that the number of cells is a function of 5 (N = n(5)) 

l=n{5)5 d . (9) 

Taking logarithms of both sides of (|D|) and solving for d, we obtain the standard definition 
of dimension: 

,= -!ff) (io) 

Suppose now that we want to measure the dimension of a Cantor set, constructed as 
follows (see figure): 

We take each segment (figure 7A), divide it into three equal parts, remove the middle 
segment and replace it with two other segments of the same length (figure 7B), repeating 
this procedure ad infinitum (figure 7C, ..), one gets (figure 7D). 

To cover the Cantor set (figure tal n) with segments it would be necessary to use an 
infinitesimal e, implying that 5 — > 0. Hence the analogue of equation ( |r0|) would be 

d = -hm ; ; " 11 

In (5) v ' 

which is the Hausdorff dimension. 

In the Cantor set above we can easily determine the value for the limit on (|TT|): 

ln(4) , N 

d = H-~1.26. 12 

ln G) 



A. 2 Fractal Dimension of Boundaries 

Consider a rf-dimensional finite region R, and a subset S C R. How can we determine 
the fractal dimension of S, given some criteria as to whether a point of R belongs to S 
or not? 

One possible way to determine the fractal dimension of S is as follows. Without loss of 
generality, let us suppose that R is a hypercube of side a. Dividing R into N hypercubic 
"cells" of side e, and evaluating the number of cells N s that contain points belonging to 
the subset S we expect that for 5 < 1 we have a good approximation to the exact result 
given in (III]). 

However the number of cells 

N = 5- d (13) 



16 



rapidly increases as 5 — > 0, and this approach becomes impractical. An alternative ap- 
proach is that of selecting N* random cells in R and counting the number of cells N* that 
have points belonging to S. We expect that 



N* ~^ N 

as N* becomes statistically sound. Supposing that, indeed, 

N* N 

we have 

ln(i\g ]n(N*/N*) ln(iV) 
In 5 In 5 In 5 



Defining 



we have 

'N* 

ln(-M =a In 5. 



(14) 

(15) 
(16) 



a = d-d f = \ (17) 



Both 5 and N* /N* can be measured, allowing us to determine a and the fractal dimension 
df of S. 

In the previous discussion, little has been said about how to determine whether or not 
a point belongs to the subset S. In our specific case of differential equations, the following 
criterion is used. For a given cell, a number of points within the cell is chosen, and the 
trajectories passing through those points are followed. If these trajectories have different 
attractors, the cell is said to belong to S. 
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Figure Captions 



Figure 1: A plot of the Lorenz system (|J) using Nsolve. 

Figure 2: \n(N/N ic ) as a function of lne. N ic = 10 4 and HT5 < e < 10~ 4 . The 
slope of the curve, 0.13 ± 0.01, is the exponent a in ([3D defining the fractal dimension: 
D = 2.87±0.01. 

Figure 3:As for Figure 2, but with 10 -8 < e < 10~ 7 . The slope of the best-fit straight 
line is a = 0.12 ± 0.01. 

Figure 4:The output from Fdimension applied to the non-chaotic (non-fractal) Lorenz 
system (a = 10, b = 8/3 and R = 10), with 5 x 10~3 < e < 9 x 10" 3 . The slope of the 
best-fit straight line is a = 1.02 ± 0.05. 

Figure 5:An example of the use of View for a non-fractal system. 100 (random) initial 
conditions have been taken in the phase space interval [x = — 0.1. .0.1, y = — 0.1. .0.1, z = 
0.9. .1.1]. 

Figure 6:The result of 100 (random) initial conditions taken in the phase space in- 
terval [x = 0.9. .1.1, y = 0.9.. 1.1, z = 21. 9. .22.1], for a chaotic regime. 
Figure 7:The construction process for the Cantor Set 
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